\documentclass[UTF8]{ctexart}
\usepackage{ctex}
\usepackage{amsmath}
\usepackage{amsthm}
\usepackage{geometry}
\geometry{left=2.5cm,right=2.5cm,top=2.5cm,bottom=2.5cm}
\usepackage{amssymb}
\usepackage{indentfirst}
\usepackage{graphicx}
\usepackage{subfigure}
\usepackage{listings}
\usepackage{xcolor}
\usepackage{float}
\usepackage{algorithm}  
\usepackage{algorithmicx}  
\usepackage{longtable}
\usepackage{fancyhdr}
\usepackage{appendix}
\usepackage{enumitem}
\usepackage{abstract}
\usepackage{multirow}
\pagestyle{fancy}
\lfoot{}%这条语句可以让页码出现在下方
\theoremstyle{plain}
\newtheorem{thm}{Theorem}[section]
\newtheorem{lem}[thm]{Lemma}
\newtheorem{prop}[thm]{Proposition}
\newtheorem{cor}[thm]{Corollary}

\theoremstyle{definition}
\newtheorem{defn}{Definition}[section]

\theoremstyle{remark}
\newtheorem*{rem}{Remark}
\newtheorem{eg}{Example}[section]
\title{Math document}
\author{Shuang Hu}
\begin{document}
\maketitle
\section{为何需要多重网格}
多重网格法是一种求解线性方程组的极为特殊的方法。

关于线性方程组，算法倒是有相当多，比如传统的高斯消去法，LU分解法（$O(N^3)$）,正定矩阵的Cholesky分解（$O(N^2)$），经典迭代法（最高可达$O(N^{1.17})$）等等。相比之下，针对模型问题的多重网格法是一种相当特殊的算法，它的时间复杂度甚至可以优化到$O(N)$，这在实际计算中，意义将是非常重大的。
\section{1D模型问题}
一维模型问题可以进行如下定义：
$$
\left\{
\begin{aligned}
-u''(x)&=f(x) \quad(x\in(0,1))\\
u&=0\quad(x\in\{0,1\})
\end{aligned}
\right.
$$
对该问题设计有限差分法格式，可以得到线性方程组$Au=f$,其中$A\in\mathbb{R}^{(n-1)\times(n-1)},n=\frac{1}{h}$：
$$
A=\frac{1}{h^2}\begin{pmatrix}
2&-1&0&\cdots&0&0\\
-1&2&-1&\cdots&0&0\\
0&-1&2&\cdots&0&0\\
\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\
0&0&0&\cdots&2&-1\\
0&0&0&\cdots&-1&2\\
\end{pmatrix}
$$
$b$则与右端项$f(x)$相关。那么，最后一步的问题自然是：如何求解方程组$Au=b$。

作为引理，先可以给出矩阵A的特征值和特征向量如下:
$$
\begin{aligned}
\lambda_{k}(A)&=\frac{4}{h^{2}}\sin^{2}(\frac{kh\pi}{2})\\
\omega_{k,j}&=\sin\frac{jk\pi}{n}
\end{aligned}
$$
\section{经典迭代格式}
可以看出，当$h$取一个很小的值，这里的矩阵A是一个大规模稀疏矩阵。在数值线性代数课程中，这种大规模的稀疏矩阵往往采用古典迭代法求解。写带权Jacobi迭代格式如下：
$$
u^{(l+1)}=T_{w}u^{(l)}+c,
$$
其中
$$
T_{w}=I-\frac{wh^{2}}{2}A.
$$
根据上面所述矩阵A的特征值，可得迭代矩阵$T_{w}$的特征值为
$$
\lambda_{k}(T_{w})=1-2w\sin^{2}\frac{k\pi}{2n},k=1,2,\cdots,n-1
$$
迭代法的收敛速度取决于迭代矩阵的谱半径。然而，可以想象，当$n$很大时，$T_{w}$的模最大特征值肯定相当接近1，这会使得迭代收敛的速度大大降低。如果要让$T_{w}$的谱半径尽可能小，或许需要$n$比较小才行，但这对方程求解的精度又是一个巨大打击...

那么，能不能在边值问题的数值求解的过程中，对步长$h$进行一些动态修正，从而兼顾迭代速度和求解精度呢？

这就是多重网格法的主要思路。
\section{格点上的Fourier分析}
观察上一节中$\lambda_{k}(T_{w})$的形式，可以发现主要的麻烦出在$k$很小的时候。那么一个很自然的想法呼之欲出：能不能把$k<\frac{n}{2}$和$k\ge\frac{n}{2}$分开考虑，尽可能避免$k$小的情形，接纳$k$大的情形？

这时自然就绕不开关于特征向量$w_{k}$的讨论。

由于$A$是一个实对称矩阵，根据线性代数知识可知$w_{k}$彼此正交，并组成一组基。

把迭代法求解的误差向量$e$按这组基进行分解，由于$\lambda_{k}$不同，可知$e$在不同的方向$w_{k}$上的衰减速度可能会有不同。如果$k>\frac{n}{2}$，我们称对应的$w_{k}$为\textbf{高频波}，而如果$k<\frac{n}{2}$，则称$w_{k}$为\textbf{低频波}。

可以预料的是，针对方程$Au=f$采用带权Jacobi迭代的过程中，高频波会迅速衰减，但低频波衰减速度则会很慢。很自然的想法是，我们在对步长进行动态修正时，能不能尽可能让原有的低频波转换为高频波？
\section{网格处理：限制和延展}
首先给出一个符号：
$$
w_{k,j}^{h}:=\sin(x_{j}k\pi)
$$
该符号表达的是网格密度取$h$时对应$w_{k}$的第$j$个分量，其中$x_{j}=\frac{j}{n},n=\frac{1}{h}$。

则有下面的引理成立：
$$
w_{k,2j}^{h}=w_{k,j}^{2h}，
$$
并且对于$k\in[\frac{n}{4},\frac{n}{2})$，网格密度为$h$时$w_{k}$为低频波，但网格密度为$2h$时即为高频波。

这就给了我们一个不错的思路：有没有可能先进行几步迭代，消去高频波，然后通过加粗网格，把低频波转化为高频波后消去，从而对解进行修正？

这么做的时候，首先涉及到的问题就是如何把细网格迭代过程中产生的残差$r=Av-f$转移到粗网格上，以及将粗网格上的残差转移到细网格上。

这两件事由下面两个算子达成：
$$
I_{h}^{2h}v^{h}=v^{2h}
$$
其中$v_{j}^{2h}=\frac{1}{4}(v_{2j-1}^{h}+2v_{2j}^{h}+v_{2j+1}^{h})$
$$
I_{2h}^{h}v^{2h}=v^{h}
$$
其中$v_{2j}^{h}=v_{j}^{2h}$,$v_{2j+1}^{h}=\frac{1}{2}(v_{j}^{2h}+v_{j+1}^{2h})$

第一个算子称为限制算子，而第二个算子则称为延展算子。
\section{多重网格算法}
有了以上的准备，就可以设计出多重网格算法。个人理解，这些算法体现了预估-校正的思路，也就是先对细网格进行几步迭代，得出求解的余项，再利用粗网格对余项进行处理，从而对所得的解进行一些校正。

算法见课本《A Multigrid Tutorial》（William L.Briggs）的page40~43
\section{程序实现时碰到的一个bug以及反思}

我在编写本次大作业的程序实现时，遇到了一个bug。由于此bug主要出现在数学推导部分，故不在设计文档中记录，而转至数学文档。

为了编写程序方便起见，我把原问题的方程组改了一改，事实上就是两边同时乘了系数$h^2$，使得系数矩阵A的元素与h无关。如下所示：

$$
A=\begin{pmatrix}
2&-1&0&\cdots&0&0\\
-1&2&-1&\cdots&0&0\\
0&-1&2&\cdots&0&0\\
\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\
0&0&0&\cdots&2&-1\\
0&0&0&\cdots&-1&2\\
\end{pmatrix}
$$
$$
b=\begin{pmatrix}
h^{2}f(x_{1})+\alpha\\
h^{2}f(x_{2})\\
\vdots\\
h^{2}f(x_{m-1})\\
h^{2}f(x_{m})+\beta\\
\end{pmatrix}
$$

我此时认为既然方程等价，算法本身自然没有修改的道理。从而我按讲义上的VCycle算法形成了程序，结果居然不收敛？

我把这个式子和讲义上的形式做了仔细对比后，发现问题出在残差上。讲义上残差$e$定义为$b-Au$，这对于不调整步长的情形，当然成立。但多重网格程序中，步长$h$会动态调整。步长为$2h$时，用讲义上的方程，对应求得的残差为$F_{2h}-A_{2h}u$，那对于我们的方程，两边乘以$h^2$，残差则为$\frac{1}{4}(\tilde{b}-A\tilde{u})$，较之真实残差刚好少了$\frac{1}{4}$!

因此，我对程序进行了修正，体现在每次计算完残差后都先乘4。再对改完的程序做测试，算法收敛了。
\end{document}